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We report a numerical analysis of Anderson localization in a model 
>■£ , of a doped semiconductor. The model incorporates the disorder arising 

from the random spatial distribution of the donor impurities and takes 
account of the electron-electron interactions between the carriers using 
density functional theory in the local density approximation. Preliminary 
results suggest that the model exhibits a metal-insulator transition. 



1 Introduction 



In semiconductors a zero temperature metal-insulator transition is observed as a 
function of doping concentration. For samples with concentrations below a crit- 

q . ical concentration, the conductivity extrapolated to zero temperature is found 

to be zero. For samples with concentrations exceeding this critical concentra- 
tion, the zero temperature limit of the conductivity is finite.^ One well studied 
pg ' example is phosphor doped silicon (Si:P) (see Ref. [2] for a recent review). The 

J> , relative importance of the roles that electron-electron interactions and disorder 

play in this transition is still not clear. The Coulomb interaction between the 
electrons leads us to expect that the impurity band is split into upper and lower 
Hubbard bands and that the transition is associated with a closing of the Hub- 
bard gap. However, this ignores the effect of the disorder that arises from the 
random spatial distribution of the donor impurities and the possibility of Ander- 
son localization. El This paper is a preliminary report of numerical simulations 
designed to address this issue. 



2 Model 



As a simple model of a doped semiconductor we consider an effective medium 
with electron effective mass m* and dielectric constant e r equal to those of the 
host semiconductor crystal. In this effective medium TV donor impurities are 
randomly distributed in space. Since we have mind phosphor in silicon, we as- 
sume that each donor supplies one electron and as a result has a net charge of 
+e. This is the only property of the donor which enters our model. There are 
an equal number of electrons so that the total charge is zero. The electrons 
interact with the donors through the Coulomb interaction. The random spa- 
tial distribution of the donors thus produces a random potential in which the 



electrons move. At the same time the electrons interact with each other via the 
Coulomb interaction. The Hamiltonian of this system is 



1 N , 1 N 1 1 1 

c 1=1 1,1=1 ' ' l=£j ' J ' 

Here, Hartree atomic units are used. The positions of the donor impurities 
are denoted by {Ri}- The first term is the kinetic energy of the electrons, the 
second term describes the interaction of the electrons with the donor impurities, 
and the third term describes the interaction between the electrons. A fourth 
term describing the mutual Coulomb interaction between the donor impurities 
should also be included if the correct total energy of the system is required. 
However, since the positions of the donor impurities, while random, are fixed, 
this contribution to the energy does not play any role in the following discussion 
and is, therefore, omitted. 

To deal with the electron-electron interaction we use density functional 
theory SI and solve the Kohn-Sham equations ^ that describe an auxiliary one- 
electron problem that has the same ground state density as the interacting 
problem of Eq. ([!]). The Kohn-Sham equations are 



(-^^ 2 + V eS [n(r)]U i (r) = e i <p i (r) (i = 1,. . .,N) , 



(2) 



where 



ff [n (r)] = ^E^ + iJ^> i^i + V XC [n (r)] . (3) 



The number density of the electrons is 

N 

n(r)=Y,\<l>i(r)\ 2 ■ (4) 

i=l 

Periodic boundary conditions are imposed. In this model the dependence on the 
medium enters only through the effective mass and dielectric constant. Having 
in mind silicon as the host semiconductor we set 

m* = 0.32 , e r = 12.0 . (5) 

The exchange-correlation potential appearing in the Kohn-Sham equations 
is given by the functional derivative of the exchange-correlation energy with 
respect to the number density of electrons 

on [r) 

While in principle the Kohn-Sham equations are exact, in practice the exact 
form of the exchange-correlation potential is not known and an approximation 
is required. In this work, we use the local density approximation (LDA) in 
which the functional is approximated as 

E XC « ££DA s I d 3 rexc (n (r)) n (r) _ (?) 



In this preliminary work, we assume complete spin polarization^ We use the 
form of exc given in Eq. (2) of Ref. [5] (with spin polarization £ = 1) though 
with the parameter values given in Ref. [7] rather than Ref. |H1 

In the literature expressions for the exchange-correlation potential are given 
for electrons in free space whereas we are considering an effective medium. To 
map the expressions in the literature to the formulae we require here, we re-scale 
lengths and energies according to the formulae 

r = (m*Je r )r, E = (s 2 r / 'm* e ) E . (8) 

After this re-scaling we have 



where 



777 _ 

Vxc[n(r)] = ^-V xc [h(f)} , (9) 

F 

~ r 



n(f) = (e r /m* e ) 3 n(r) , (10) 



and Vxc is the exchange-correlation potential found in the literature. 

In Hartree atomic units the unit of length is the Bohr a,Q = e^h? / [irm e e 2 ) sj 
5.292 x 10 _11 m. We use this as the unit of length throughout. The simulations 
are performed by generating an ensemble of cubic samples of linear dimension 
L. In each sample, donor impurities are randomly distributed on a cubic lattice 
with lattice constant 36 Bohr. This prevents impurities being positioned un- 
physically close together by chance. (We do not allow two impurities to be at 
the same site on this lattice.) The volume of the system is V = L 3 and the 
donor concentration njj = N/V. 

For numerical purposes the continuous description above is replaced by a 
discrete description on a real-space grid with spacing a. Derivatives are re- 
placed by next nearest neighbor finite difference approximations. The resulting 
matrices and vectors have dimension equal to the number of grid points (L/a) . 
When calculating the results shown below we used a grid spacing of 

a =18 Bohr. (11) 

For comparison the effective Bohr radius for an electron in the conduction band 
of silicon is aj = £ r /m e Bohr rj 37.5 Bohr. The potential due to the positive 
donor impurities is calculated by expressing the charge density of the impurities 
as a Fourier series. A cut-off is imposed on the wavenumbers so that the number 
of terms in this series is equal to the number of points of the real-space grid. In 
effect, this replaces the delta- functions of the charge density of the impurities 
with an approximate smooth charge density. Poisson's equation is solved exactly 
for this approximate density and the corresponding potential obtained using an 
inverse Fourier transform. This calculation need only be performed once for a 
given impurity configuration. The Hartree like term in Eq. ([3]) is evaluated in 
a similar way. This latter calculation needs to be repeated for each iteration. 
Within the LDA the real space finite difference approximation of the Kohn- 
Sham equations yields a Hamiltonian matrix that is sparse. The N Kohn-Sham 
orbitals of lowest energy are found using the JADAMILU sparse matrix library. El 
The self consistent solution of the Kohn-Sham equations is found by iteration. 



1 A calculation that assumes either zero or complete spin polarization is numerically less 
demanding than a calculation using the local spin density approximation but excludes the 
possibility of a Hubbard gap from consideration. Wc do not assume zero polarization because 
double occupation of the impurities is unlikely in the localized phase. 
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Figure 1: The value of generalized multi- fractal exponent ao as a function of donor impurity 
concentration for L = 432, 540, and 648 Bohr. The number of samples simulated varies 
between 494 and 1004. 



3 Result 

To determine the nature, localized or extended, of the Kohn-Sham wavefunc- 
tions we use multi-fractal finite size scaling ED LBS (MFSS) of the wavefunction 
intensity distribution. The system is divided into boxes of linear dimension I 
and the Kohn-Sham wavefunction intensity is coarse grained 
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' kth box 

From these coarse grained intensities a random variable a is defined 

hlyU 

a =h7A' 

where A is the ratio of box and system sizes 

>4- 



(12) 



(13) 



(14) 



As described in Refs. and [TOJ for Anderson's model of localization^! the 
distribution of a is scale invariant at the Anderson transition provided A is held 
fixed. The distribution shifts to smaller (larger) values in the metallic (localized) 
phases as the systems size increases. We expect similar considerations to apply 
to the intensities of the Kohn-Sham wavefunctions. In what follows we focus on 
the generalized multi-fractal exponent ao for the highest occupied Kohn-Sham 
orbital. We refer the reader to Eqs. (6), (7) and (19) of Ref.lTOlfor the definition 
of this quantity. The ensemble average in the definition given there is estimated 
using an average over samples in the usual way. The precision of this estimate 
is determined using the formulae given in Table II of Ref. [TO] 



In Fig [T] we plot the donor concentration dependence of So for three different 
system sizes. Least squares fits of third order polynomials have been made 
to the data for each system size. For low concentration behavior typical of 
Anderson localized wavefunctions is seen, i.e. we see shifts to larger values as the 
system size increases. For high concentrations the system size dependence is less 
pronounced but there seems to be a shift to smaller values with increasing system 
size. (Note there is a lower bound of So > 3 set by wavefunction normalization.) 
The polynomial fits show a common crossing point around 1.4 x 10 _7 Bohr _ . 
The data, though certainly very preliminary, suggest a metal-insulator transition 
near this value. For comparison the experimental critical concentration for 
the metal-insulator transition in Si:P is ~ 5.2 x 10~ 7 Bohr~ .^ It is possible 
that this discrepancy arises from our simplifying assumption of complete spin 
polarization. 

Examination of the electron density (not shown here) for the highest oc- 
cupied Kohn-Sham wavefunction of typical samples reveals the following. For 
low concentrations ~ 0.6 x 10 _7 Bohr~ the Kohn-Sham wavefunctions resemble 
molecular orbitals on clusters of two or three impurities. As the concentra- 
tion is increased towards ~ 1.4 x lCP 7 Bohr~ 3 the Kohn-Sham wavefunctions 
spread out over more impurities but remain localized. For higher concentra- 
tion ~ 1.8x 10~ 7 Bohr~ the Kohn-Sham wavefunctions are extended across the 
entire sample. 

4 Conclusion 

We have presented a model for a doped semiconductor that includes both the 
random spatial distribution of the donor impurities and the interaction between 
the carrier electrons, and whose numerical analysis allows the metal-insulator 
transition to be studied. Preliminary results indicate that the model does ex- 
hibit a metal-insulator transition. Since in Eq. (TTJ the properties of the semi- 
conductor enter only through the effective mass and relative dielectric constant, 
we automatically recover the behavior described in Fig. 1 of Edwards and 
SienkoEEQ that to a good approximation the critical concentrations in various 
semiconductors obey 

n c a o = constant . (15) 

However, our preliminary estimate of this constant, rj 0.19, differs significantly 
from the observed value of ~ 0.26. As mentioned already, it is possible that this 
discrepancy will be resolved when our calculations are extended to properly take 
account of the electron spin. 

In future work we hope to shed light on the relative importance of the roles 
of disorder and interaction in this transition by making a careful comparison of 
the critical properties of this model with the known properties of the Anderson 
transition in Anderson's model of localization for non-interacting electrons. Due 
to space limitations we also defer a discussion of previous workE^l E3 EH. 
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